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Preliminary large-eddy simulations of flow around 
a NACA 4412 airfoil using unstructured grids 

By Kenneth Jansen 


1. Motivation and objectives 

Large-eddy simulation (LES) has matured to the point where application to com- 
plex flows is desirable. The extension to higher Reynolds numbers leads to an im- 
practical number of grid points with existing structured-grid methods. Furthermore, 
most real world flows are rather difficult to represent geometrically with structured 
grids. Unstructured-grid methods offer a release from both of these constraints. 
However, just as it took many years for structured-grid methods to be well un- 
derstood and reliable tools for LES, unstructured-grid methods must be carefully 
studied before we can expect them to attain their full potential. 

In the past two years, important building blocks have been put into place mak- 
ing possible a careful study of LES on unstructured grids. The first building block 
was an efficient mesh generator which allowed the placement of points according to 
smooth variation of physical length scales. This variation of length scales is in all 
three directions independently, which allows a large reduction in points when com- 
pared to structured-grid methods, which can only vary length scales in one direction 
at a time. The second building block was the development of a dynamic model ap- 
propriate for unstructured grids. The principle obstacle was the development of 
an unstructured-grid filtering operator. New filtering operators were developed in 
Jansen (1994). In the past year, some of these filters have been implemented into 
a highly parallelized finite element code based on the Galerkin/least-squares finite 
element method (see Jansen et al. (1993) and Johan et al. (1992)). 

We have chosen the NACA 4412 airfoil at maximum lift as the first simulation 
for a variety of reasons. First, it is a problem of significant interest since it would 
be the first LES of an aircraft component. Second, this flow has been the subject 
of three experimental studies (Coles and Wadcock (1979), Hasting and Williams 
(1987), and Wadcock (1987)). The first study found the maximum lift angle to be 
13.87 . The later studies found the angle to be 12°. Wadcock reports in the later 
study that the early data agree very well with his new data at 12°, suggesting that 
the early experiment suffered from a non-parallel mean flow in the Caltech wind 
tunnel. It should be pointed out that the Reynolds- averaged simulations are usually 
run at 13.87° and do not agree with the data when run at 12° as will be shown later 
in this study. It is hoped that LES can clarify this controversy. The third reason for 
considering this flow is the variety of flow features which provide an important test of 
the dynamic model. Starting from the nose where the flow stagnates, thin laminar 
boundary layers are formed in a very favorable pressure gradient. This pressure 
gradient soon turns adverse, driving the flow toward a leading edge separation. 
Only the onset of turbulence can cause the flow to remain attached or to reattach 
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if it did separate. The persistent adverse pressure gradient eventually drives the 
turbulent flow to separate in the last 20 percent of chord. The separation bubble is 
closed near the trailing edge as the retarded upper surface boundary layer interacts 
with the very thin lower surface boundary layer. The large difference in boundary 
layers creates a challenging wake to simulate. Only the dynamic model can be 
expected to perform satisfactorily in this variety of situations: from the laminar 
regions where it must not modify the flow at all to the turbulent boundary layers 
and wake where it must represent a wide variety of subgrid-scale structures. 

The flow configuration we have chosen is that of Wadcock (1987) at Reynolds 
number based on chord Re c = u^c/u = 1.64 x 10 6 , Mach number M = 0.2, and 
12° angle of attack. 

2. Accomplishments 

2.1 Dynamic model implemented and tested 

The only obstacle to implementing a dynamic model on unstructured grids is 
extension of the filtering operator. Four filtering operators were proposed in Jansen 
(1994). Two of these models were implemented and compared using a simple ana- 
lytic velocity field for which the filtered values can be determined exactly. From this 
test, the generalized top-hat was found to be the most accurate, and all subsequent 
calculations have been carried out using this filter. 

2.2 Simulations 

A series of simulations has been performed in the last year to develop experi- 
ence with this new approach. The first simulation was intentionally very coarse 
as we hoped to improve the mesh selectively and develop an understanding of the 
sensitivity of the solution to the grid improvements. 

2.2.1 First simulation 

The first simulation was performed on a very coarse mesh. The near- wall grid 
did not attempt to resolve the near-wall layer accurately in the first 20 percent of 
chord and only marginally resolved the remaining flow (A* = 300, A* = 80) at 
the wall. The grid was coarsened in the streamwise and spanwise directions coming 
off the wall as suggested by Chapman (1979). The resulting coefficient of pressure 
distribution was reasonably well predicted on this mesh (see Fig. 1), but the velocity 
profiles showed poor agreement with the experiment. 

2.2.2 Improvement of outer layer 

Careful scrutiny of the mesh revealed that the strategy of coarsening in the 
streamwise direction coming off the wall was inappropriate at this Reynolds number. 
The inner-layer spacing, A, scales on wall units. 

A v y/2 

c u r c Re cy /Cf 

Here, i/, is the kinematic viscosity, c is the chord length, and u r is the friction 
velocity defined to be the square root of the coefficient of friction, C/, over two. 
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Figure 1 . Coefficient of pressure along the airfoil surface. LES , Hastings 

and Williams o , Wadcock (1987) □ . 

The outer-layer spacing scales on the boundary layer thickness, S 99 . It is reasonable 
to expect the large eddies in the outer part of the boundary layer to be of order S 99 , 
and therefore the outer-layer spacing, in all directions, should never exceed 



By using Wadcock s experimental data for the C/ and S 99 , one can compare these 
two resolution restrictions as is done in Fig. 2. This figure contains three curves. 
The solid curve describes the variation of a 200 wall-unit spacing (which can be 
associated with the streamwise spacing near the wall) over the upper surface where 
the boundary layer is attached. The dashed curve describes the same variation 
of 50 wall units (which can be associated with spanwise spacing near the wall). 
The chain dash curve is the outer-layer spacing as described above. Several points 
can be made in this figure. First, all three curves change by over an order of 
magnitude from the tip to the tail region. This illustrates how an unstructured grid 
saves points by matching resolution to the local changes in the length scales in the 
streamwise direction. For example, a structured grid would be forced to carry the 
fine spanwise resolution required near the nose through the entire domain. Second, 
when comparing the near-wall spanwise resolution to the outer-layer resolution^ 
it is clear that coarsening the spanwise resolution as the distance from the wall 
increases is justified. The final point, apparent from this figure, is that coarsening 
of the streamwise resolution in the outer layer is not justified. In fact, over much 
of the airfoil surface the outer-layer grid resolution is more restrictive than the 
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FIGURE 2. Comparison of length scales over the airfoil surface. 200 wall units 

(streamwise near- wall spacing) , Sgg / 5 (outer-layer spacing) , 50 wall 

units (spanwise near-wall spacing) . 

inner-layer resolution. The choices of 200 wall units and 5 points per boundary 
layer thickness are somewhat arbitrary, but they axe believed to be comparable in 
their degree of coarseness. It is interesting to observe that the crossover between 
these two curves corresponds to Re r = (u r 6gg / u) = 1000. Therefore, when above 
1000, the inner-layer resolution is the most restrictive. Otherwise, the outer-layer 
resolution is the most restrictive. Only at higher Reynolds numbers will coarsening 
in the streamwise direction be justified. 

Considering the above discussion, a new mesh was made where the coarsening 
of the streamwise spacing was delayed until outside of the boundary layer. This 
resulted in a mesh with nearly twice as many points as the previous simulation. It 
also resulted in a rather dramatic change in the early boundary layer structure. It 
seems that the improved resolution of the outer layer allowed a better resolution of 
the leading edge separation. The new simulation led to a train of spanwise coherent 
vortices. These vortices broke down into turbulence at about 10 percent of chord. 

The persistence of the spanwise coherent vortices was not in line with the exper- 
iments which were all tripped. Some evidence as to the importance of the tripping 
can be seen in Fig. 3 where we compare the surface coefficient of pressure distribu- 
tion from the free transition simulation to two experimental data sets from Hastings 
and Williams (1987). The square data set was taken without a transition strip while 
the circle data set was taken with a transition strip. Our simulation shows rather 
good agreement with the free transition. Unfortunately, all velocity and Reynolds 
stress data were taken with the transition strip in place and agreement with these 
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Figure 3. Coefficient of pressure along the airfoil surface. LES , Hastings 

and Williams without trip □ , Hastings and Williams tripped o . 


quantities is substantially worse than with C p . 

2.2.S Grid refinement study of the nose 

The dramatic change of the flow with changed resolution indicated a need for 
further refinement in the nose region. At the same time we also hoped to model the 
transition through a steady blowing pattern as shown in Fig. 4. A shape that could 
be easily resolved was chosen. Therefore, we could be certain that any sensitivity 
to grid refinement would be associated with the turbulence structures responding 
to the blowing and not the resolution of the blowing itself. 

A new mesh was generated where the streamwise and spanwise resolution were 
improved by a factor of two everywhere on the upper surface. The normal spacing 
was improved at the wall by a factor of two as well, but this did not lead to a 
doubling of points in this direction due to the stretching. The spanwise domain 
was cut in half (from 0.05c to 0.025c) for this simulation. Therefore, the number of 
points approximately doubled rather than a quadrupling. 

There was again a rather dramatic change in the solution and so another mesh 
was generated. This mesh again improved the streamwise and spanwise resolution 
by a factor of two, although, this time, only in the first 5 percent of chord. The three 
surface meshes of the first 10 percent of chord are shown in Fig. 5. The velocity 
profiles in the first 5 percent of chord are shown in Fig. 6. For this forcing pattern, 
the flow is nearly spanwise- and streamwise-resolution independent. 
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FIGURE 4. Elevation plot of the steady jet normal to the airfoil surface. The 
actual grid is shown to confirm resolution. 




FIGURE 5. Surface meshes near the leading edge (0.0 < xj c < 0.1). Mesh (b) has 
been refined by a factor of 2, both spanwise and streamwise, from mesh (a). The 
spanwise domain is also halved. Mesh (c) has been refined by a factor of 2, both 
spanwise and streamwise, from mesh (b) in the first 5 percent of chord. 
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Figure 6. Profiles of tangential velocity component at various positions along 
the airfoil surface (x/c = 0.005,0.01,0.02,0.03,0.04,0.05). Solutions correspond to 
grids from previous figure; mesh (a) , mesh (b) , mesh (c) 



Figure 7. 


Boundary layer parameters (lower curves and data are momentum 


thickness ( 6 ), higher curves and data are displacement thickness (6*)) along the 

airfoil surface. Solutions correspond to meshes from Fig. 5; mesh (a) , mesh 

(b) , mesh (c) , Wadcock (1987) (o ). 
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2.3 More accurate transition 

While it was useful to obtain a grid-independent solution at the forcing prescribed, 
the final solution does not agree with experiment, as can be seen in Fig. 7 . Here, 
the momentum and displacement thickness of the grid-independent calculations 
can be seen to be substantially greater than the experiments at the first available 
datum point. The discrepancy seems to be associated with the laminar separation 
at 1 percent of chord. The simulation suggests a transition in the free shear layer, 
followed by a turbulent reattachment. This mode of transition seems to give the 
flow a large jump in momentum and displacement thickness. The experiment did 
not seem take this route to transition. For this reason a more careful study of 
transition is currently underway. 

Wadcock used a strip of tape with serrations cut into the edge on the upstream 
side. The serrated tape can be modeled in a coarse fashion by our current simulation 
as can be seen in Fig. 8. The tape is effectively a forward facing step (with serrations) 
of height 699/4, followed by a backward facing step. Calculations are underway with 
this modification. 



FIGURE 8. A transition strip is modeled geometrically by applying a no slip 
boundary condition to the nodes which form a surface of height, shape, and position 
equivalent to Wadcock’s serrated tape which was applied to the airfoil surface. 

2.4 Reynolds- averaged simulations 

Reynolds- averaged Navier-Stokes simulations (RANS) have not shown good agree- 
ment with the experimental data. However, given the cost of LES, they can be a 
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helpful tool for suggesting sensitivity to changes of basic flow parameters since they 
require so little computational time. While the results are not expected to be quan- 
titatively correct, trends can at least be suggested by RANS and later confirmed 
by LES. 

A series of RANS calculations was performed to chart various trends in this 
flow. The RANS calculations used the commonly accepted NASA code (INS2D) of 
Rogers (1991) and employed a k-u model from Menter (1994). First, the effect 
of angle-of-attack and wind tunnel walls are compared in Figures 9 and 10. The 
boundary condition on the wind tunnel walls is a slip condition. This accounts for 
the blockage of the walls without requiring resolution of the boundary layers on 
them. The effects are compared together because it is common among the RANS 
modeling community to adjust the angle-of-attack of free air calculations to account 
for the walls. Figure 10 suggests that the flattening of the C p near the trailing edge 
(which is associated with the large separation there) is affected strongly by angle- 
of-attack and only weakly by the wind tunnel walls. The 13.87° angle-of-attack 
cannot be justified with the hope of accounting for the effects of the wind tunnel 
walls in free air calculations. 

The second trend studied with the RANS code was the effect of transition posi- 
tion. When the RANS code was run with the transition point fixed at the position 
of Wadcock’s strip, a leading edge separation developed on sufficiently fine meshes. 
Once beyond the transition point, the flow reattaches. This provides an independent 
verification of the results observed in the LES. 

3. Future plans 

3-1 Grid- independent solution of flow with a transition strip 

The calculation using the transition strip described above will be continued and 
checked for grid dependence. It should be noted that grid independence can only be 
achieved beyond a short distance downstream of the transition strip. True grid inde- 
pendence of the strip and transition itself is probably too expensive to be practical, 
even with an unstructured grid. It may be necessary to provide small disturbances 
upstream of the strip to mimic the interaction of freestream turbulence with the 
strip. This capability has been implemented and tested in the code using a wall jet 
with spatial and temporal variation. 

3.2 Inclusion of the wind tunnel walls 

The RANS studies indicated a moderate effect of the wind tunnel walls on the 
solution. Future simulations will be done with a slip boundary conditions on the 
wind tunnel walls. Meshes have already been generated for this purpose as can be 
seen in Fig. 11. 

3.8 Higher order methods 

Given the number of points that are required to get a grid-independent solution, 
it seems clear that higher order methods should be explored. This is relatively easy, 
but non- trivial, to do with the finite element method. There are two benefits to 
higher order methods besides the obvious one of higher accuracy. First, the higher 



FIGURE 9. Coefficient of pressure along the airfoil surface from RANS simulations. 

12° in free air , 12° with walls , 13.87° in free air , 13.87° with 

walls , Hastings and Williams (1987) o , Wadcock (1987) □ . 



FIGURE 10. Closeup of the trailing-edge region (same legend as above). 
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Figure 11. The crosssectional plane of an unstructured mesh which accounts for 
the inviscid effects of the wind tunnel walls. 

order methods will have a more complete representation of the residual error of the 
discrete approximation and, therefore, the scheme will be less dissipative. Second, 
alternative filters, described in Jansen (1994), can be implemented and tested. It 
is difficult to predict at this time if the method will lose computational efficiency 
when extended to higher order. 

3.4 Expanded spanwise domain 

Once we are satisfied with the solution in the region of the nose, we will have to 
consider carefully the effect of the narrow spanwise domain on our solution. It is 
likely that as we predict a larger separation at the trailing edge, the effect of the 
narrow domain will become more acute. Strategies are being developed to expand 
only the portion of the domain suffering from a narrow box. If these strategies work 
a large number of points can be saved. 
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